{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![alt text](https://github.com/callysto/callysto-sample-notebooks/blob/master/notebooks/images/Callysto_Notebook-Banner_Top_06.06.18.jpg?raw=true)\n",
    "\n",
    "# Curve Fitting in Python \n",
    "\n",
    "Curve fitting is the buisiness of slapping a line (or curve) which best passes through data. By fitting functions through data, we are able to use these fitted functions in order to extrapolate and interpolate new values that we may not have measured. Let's take a look at how we can do that in Python. To begin, we import all the functions we shall need. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numerical python package to allow us to do math quickly \n",
    "import numpy as np\n",
    "# Ploting library \n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "# SciPy - package which allows us to fit curves to data\n",
    "from scipy.optimize import curve_fit\n",
    "# Helper function to make graphs without clutterin up this notebook\n",
    "from plotter import plot_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin our function fitting - we first need to have some data in order to fit! As this is simply a tutorial, let's just generate some data in order to test our fitting functions. This way we'll _know_ exactly what parameters our curve fitting functions should find. In this case, we define a linear equation of form\n",
    "\n",
    "$$ f(x) = m \\; x + b $$\n",
    "\n",
    "in Python below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_function(x, m, b):\n",
    "    return m * x + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fantastic! Now, let's generate some $x$ points for our data. To do this we will use the `np.linspace` function which creates an equaly spaced set of data. We do this below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here, we're creating 15 equally spaced numbers from 0 to 10\n",
    "x_data = np.linspace(0,10,15)\n",
    "print(x_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our $x$ values, let's _generate_ a set of $y$ values using our `linear_function` Python function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_data = linear_function(x =x_data, m = 0.5, b=-1 )\n",
    "print(y_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where we can test to see if this worked by plotting our data below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(x_data, y_data)\n",
    "plt.xlabel(\"$x$\", size =15)\n",
    "plt.ylabel(\"$y$\", size = 15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where that's fantastic! We've generated a set of linear points. However, thos points are perfectly in line! We'll certainly never see that in an actual measurement scenario. Let's add some noise to our $y$ points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we're adding normally distrubuted noise to our distributions\n",
    "noise_strength = 0.1\n",
    "y_noise = noise_strength*np.random.normal(size=x_data.size)\n",
    "y_data = y_data + y_noise\n",
    "plt.scatter(x_data, y_data)\n",
    "plt.xlabel(\"$x$\", size =15)\n",
    "plt.ylabel(\"$y$\", size = 15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There we go! Now that's a little more realistic for function fitting. Now that we've generated some data, let's see how we fit a function using Python.\n",
    "\n",
    "## The `curve_fit` Function\n",
    "We will be using the `curve_fit` function that we imported earlier. Curve fit takes the following data:\n",
    "1. A Python function which describes the function we'd like to fit, in this example `linear_function(x,m,b)`\n",
    "2. A set of $x$ data points\n",
    "3. A set of $y$ data points \n",
    "    * We need to have the same number of points in $x$ and $y$\n",
    "    \n",
    "Once we call this function, `curve_fit` will return two things:\n",
    "1. `values`: A list of best fit parameters, in our case, it will return a list of `[m, b]`\n",
    "2. `fit_quality`: An estimate of the variance associated with our fitted function, and the data we used to generate it. These values will give us uncertainty surrounding our best fit parameters\n",
    "\n",
    "Let's take a look at that function in action\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we're fitting a curve to linear_function using x_data and y_data\n",
    "values, fit_quality  = curve_fit(linear_function, x_data, y_data)\n",
    "fit_quality = np.sqrt(np.diag(fit_quality))\n",
    "\n",
    "print(\"Our slope and intercept are:\", values)\n",
    "print(\"With uncertainty\", fit_quality)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where that's pretty good! With an original slope of $m = 0.5$, our fitted slope $m^\\prime$ came out to be $m^\\prime = 0.5 \\pm 0.05$, and our original intercept of $b=-1$ was fitted to be $b^\\prime = -1.08 \\pm 0.09$ (Note: Your values may be  different because the noise added is different every time!). Let's take a look at what the plot looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "plot_fit(linear_function, values, fit_quality, \"$f(x) = mx + b$\",x_data,y_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's all thre is to it! Play around with the amount of noise you add to the function before fitting it. How does that affect your estimates for best fit parameters as well as your estimates for uncertainty? \n",
    "\n",
    "# Fitting Non-linear functions\n",
    "\n",
    "Fitting non-linear functions to data with python is exactly the same as linear functions! Let's take a look at how we can fit data to a parabolic function of the form \n",
    "\n",
    "$$ f(x) = a\\;(x + b)^2 + c $$ \n",
    "\n",
    "in the cell below. Note we'll move a little faster this time and add our noise in the same plot!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadratic(x, a, b, c):\n",
    "    return a * (x + b) **2 + c\n",
    "\n",
    "x_quad = np.linspace(-20,20, 40)\n",
    "y_quad = quadratic(x_quad, .5, 5, 1)\n",
    "\n",
    "noise_strength = 5\n",
    "y_noise = noise_strength * np.random.normal(size=x_quad.size)\n",
    "\n",
    "y_quad = y_quad + y_noise\n",
    "\n",
    "plt.scatter(x_quad, y_quad)\n",
    "plt.xlabel(\"$x$\", size =15)\n",
    "plt.ylabel(\"$y$\", size = 15)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we fit it exaclly as we have before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values_q, fit_quality_q = curve_fit(quadratic, x_quad, y_quad)\n",
    "fit_quality_q = np.sqrt(np.diag(fit_quality_q))\n",
    "\n",
    "print(\"Values for a, b and c:\", values_q)\n",
    "print(\"Uncertainty for a, b and c:\", fit_quality_q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_fit(quadratic, values_q, fit_quality_q, \"$f(x) = a (x + b)^2 + c$\", x_quad, y_quad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How well did that fit work? Do the values we recovered from the curve fitting parameters line up with the parameters we used to generate the data set? \n",
    "\n",
    "\n",
    "## Fitting More Non-Linear Functions\n",
    "\n",
    "We can also fit even _more_ non linear functions. For example, let's generate and then fit some data generated by a normal distribution defined by\n",
    "\n",
    "$$ f(x) = \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} \\; \\exp \\left({-\\frac{(x -\\mu)^2}{2 \\sigma^2} }\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normal_function(x, sigma, mu):\n",
    "    return 1/(np.sqrt(2.0 * np.pi * sigma **2 ))* np.exp(- (x - mu)**2/(2.0 * sigma**2))\n",
    "\n",
    "x = np.linspace(-5,5,50)\n",
    "y = normal_function(x, 1,0)\n",
    "y_noise =  .01 *np.random.normal(size=x.size)\n",
    "y = y + y_noise\n",
    "plt.plot(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values, fit_quality  = curve_fit(normal_function, x, y)\n",
    "fit_quality = np.sqrt(np.diag(fit_quality))\n",
    "print(\"Sigma and mu values:\", values)\n",
    "print(\"Uncertainty in sigma and mu:\", fit_quality)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_fit(normal_function, values, fit_quality, \"Normal Distrubtion\",x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Play around with different amounts of noise and different functions to get a feel for how it works. \n",
    "\n",
    "![alt text](https://github.com/callysto/callysto-sample-notebooks/blob/master/notebooks/images/Callysto_Notebook-Banners_Bottom_06.06.18.jpg?raw=true)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "nbTranslate": {
   "displayLangs": [
    "*"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
